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We present a scalable scheme to design optimized soft pulses and pulse sequences for coherent 
control of interacting quantum many-body systems. The scheme is based on the cluster expansion 
and the time dependent perturbation theory implemented numerically. This approach offers a 
dramatic advantage in numerical efficiency, and it is also more convenient than the commonly used 
Magnus expansion, especially when dealing with higher order terms. We illustrate the scheme by 
designing 2nd-order self-refocusing 7r-pulses and a 6th-order 8-pulse refocusing sequence for a chain 
of qubits with nearest-neighbor couplings. We also discuss the performance of soft-pulse refocusing 
sequences in suppressing decoherence due to low- frequency environment. 



A control of coherent evolution of quantum systems 
is increasingly important in a number of research fields 
and applications. Such control has long been a staple in 
nuclear magnetic resonance (NMR) spectroscopy, where 
determination of the structure of complex molecules like 
proteins require the application of long sequences of pre- 
cisely designed radio- frequency (r.f.) pulses 1 . Recently, 
coherent control (CC) has emerged as an important part 
of quantum information processing (QIP), spurring nu- 
merous studies on general properties and specific design 
of pulses and pulse sequences for application in NMR- 
based 2 and other potential implementations 3 of quan- 
tum computers (QCs). This progress is closely followed 
by applications of coherent control in atomic physics 4 , 
quantum optics 5 , laser induced chemical reactions 6 , data 
communications 7 , and biomedical applications 8 . 

The precision required for QIP in particular, and for 
contemporary applications of CC in general, is achieved 
most readily using shaped (also, soft), typically narrow- 
band, pulses. When properly constructed, such pulses al- 
low excitation to be limited to a particular set of modes 
which results in better control fidelity and reduced in- 
coherent losses (e.g., heating). The latter is especially 
important for putative solid-state QC implementations 
which are proposed to operate at cryogenic temperatures. 
Additionally, as we also discuss in this work, refocusing 
with carefully designed high-order pulses and pulse se- 
quences can offer significantly better protection against 
non-resonant decoherence sources (e.g., low-frequency 
phonons) as compared to lower-order sequences. 

Over some forty years the shaped pulses were utilized 
in NMR, a number of analytical and numerical schemes 
were suggested for their design 1 . Most (although not all, 
see Ref. 9) rely on the average Hamiltonian theory, a 
perturbative scheme based on the cumulant (Magnus) 
expansion for the evolution operator. The expansion 
is done around the evolution in the applied controlling 
fields, while the chemical shifts 10 (resonant-frequency off- 
sets) and, ideally, inter-spin couplings are treated pertur- 
batively. The main drawback of the Magnus expansion 
[see Eq. (9) below] for numerics is the multiple integration 
appearing in higher orders; its use in actual calculations 
was almost always limited to quadratic order. 

The alternative scheme found in the literature is a 



simulation involving the full Hamiltonian of a quantum 
system 11 . In particular, such a calculation was done to 
optimize pulse shapes for a three-level Hamiltonian of a 
current-biased superconducting qubit 12 . This approach 
guarantees precision of custom-designed shapes, but it 
obviously lacks scalability, as the computational difficulty 
grows exponentially with the size of the system. 

In this paper, we present an efficient scalable scheme 
to design high-order soft pulses and soft-pulse sequences 
for controlling quantum many body systems. Instead 
of using the Magnus expansion or other form of the ef- 
fective Hamiltonian theory, we rely directly on the time- 
dependent perturbation theory implemented numerically. 
This allows an easy extension to higher orders (up to 9 th 
in this work), at the same time preserving the benefits 
of the cluster theorem 13 which limits the size of the sys- 
tem to be analyzed. The high order calculation allows a 
straightforward classification of pulse sequences by order 
K, the number of terms in the time dependent perturba- 
tion series for which the control remains perfect. 

To be specific, and also as an illustration, we consider 
a quantum spin chain with each spin (qubit) individu- 
ally controlled and "always-on" nearest-neighbor (n.n.) 
interactions. Assuming the Ising coupling is dominant, 
Jz J±, we construct a family of one-dimensional 7r- 
pulses 14 with different degrees of self-refocusing with re- 
spect to the J z coupling. The duration of a pulse, r, 
is chosen to be fixed, so as to allow parallel execution of 
quantum gates in different parts of the system. To reduce 
the spectral width of an arbitrary sequence of such pulses 
we also require a number of derivatives of the controlling 
fields to vanish at the ends of the cycle. 

We show that thus designed pulses work as drop-in re- 
placement of hard pulses, and compare their performance 
in several pulse sequences and composite pulses with that 
of two commonly used shapes. In particular, we present 
an eight-pulse refocusing sequence of order K — 6 for the 
quantum Ising model [refocusing errors scale as (J z r) 6 ], 
order K — 2 for general xxz model, where, in addition, 
each spin has an order K > 2 protection against phase 
decoherence due to low- frequency environment (TAB. I). 

We consider the following simplified Hamiltonian 

H{t) = H c (t) + H s + Hv(t) + H a , (1) 
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1* 1 1 
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1 2 


2* 


2* 2 3 



TABLE I: Order K determining the scaling (oc t k ) of the 
gate errors with the pulse duration r for different refocusing 
sequences. 1: a single 7r-pulse along the rc-axis on odd sites, 
"X?"; 2: two 7r-pulses along the :r-axis, odd sites only, X?X?; 
4: XiY2X!Y 2 (bars for negative pulses, subscripts denote odd 
or even sites); 8: XfYiX^YaYaX^YfX?. Asterisks mark odd- 
site refocusing. See text for description of the three models. 

with the first (main) term due to individual control fields, 

Hc(t) = \z n [v:(t) < + vy(t) *y] , (2) 

where jj, — x,y,z, are the usual Pauli matrices for 
the n-th qubit (spin) of the ID chain. The other terms 
describe the interactions between the qubits (n.n. xxz), 

H S = \ll i J n,n>«' + Jn,n>{«> + «')] , (3) 
(n,n') 

and the coupling with the oscillator thermal bath, 

tfv(t) = £ n/ ^W(t), H a = l -Y, n B%<. (4) 

In Eq. (4), A£ = A^{pi,qi) account for the possibility 
of a direct coupling of the controlling fields V£ with 
the bath variables qi, pi, while B£ = B£(pi,qi) describe 
the usual coupling of the spins with the oscillator bath. 
Already in the linear response approximation, the bath 
couplings (4) produce a frequency-dependent renormal- 
ization of the control Hamiltonian Hc(t) [Eq. (2)], as 
well as the thermal bath heating via the dissipative part 
of the corresponding response function. Both effects be- 
come more of a problem with increased spectral width of 
the controlling signals In this work we do not spec- 
ify the explicit form of the coupling Hy(t). Instead, we 
minimize the spectral width of the constructed pulses. 

Closed system: In a qubit-only system with the Ha- 
miltonian H(t) = Hc(t) + Hs, the effect of the applied 
fields is fully described by the evolution operator U(t), 

U(t) = -i[H c (t)+H s ]U(t), U(0) = L (5) 

As usual, the time-dependent perturbation theory is in- 
troduced by separating out the bare evolution operator, 

U(t) = U Q (t) R(t), U Q (t) = -iH c (t) U (t). (6) 
Then, the operator R(t) obeys the equation 

R(t) = -iH s (t)R(t), H s (t) = U^t) H s U (t), (7) 

which can be iterated to construct the standard expan- 
sion R(t) = I + -Ri(t) + R 2 (t) + . . . in powers of (tH s ), 

R n +i(t) = -iH s (t)R n (t), Ro(t) = I. (8) 



For a finite system of n qubits and a given maximum or- 
der K of the expansion, Eqs. (6), (7), and (8) are a set of 
coupled first order ordinary differential equations for the 
2" x 2™ matrices Uq, Ri, R2, ■ ■ ■ , Rk, an d can be inte- 
grated efficiently using any of the available extrapolation 
schemes. Obviously, for a given system, solving the full 
equations (5) is simpler by a factor of at least (K + 1). 
However, it is the analysis of the perturbative expansion 
that is the key for achieving the scalability of the results. 

The standard Magnus expansion can be readily ob- 
tained by integrating Eqs. (8) formally and rewriting the 
result in terms of cumulants, 

R(t) = exp(Ci + C 2 + •••), C x = -i{ d^Hsih), 

Jo 

Ci = ~j dt 2 J ' dt^Hsih^Hsfa)], ••• (9) 

Generally, the term Ck contains a fc-fold integration 
of the commutators of the rotating-frame Hamiltonian 
Hs(ti) at different time moments ti and has an order 
(tHs) h - The advantage of the cumulant expansion is that 
it does not contain the disconnected terms arising from 
different parts of the system. For an arbitrary lattice 
model of the form (3), with bonds representing the qubit 
interactions, the terms contributing to fc-th order can be 
represented graphically as connected clusters involving 
up to k lattice bonds; generally such clusters cannot have 
more than n = k + 1 vertices. Thus, to obtain the exact 
form of the expansion up to and including K-th order, 
one needs to analyze all distinct clusters with up to K + 1 
vertices. For an infinite chain with n.n. couplings, these 
are finite chains with up to K bonds and K + 1 vertices. 

The discussed cluster theorem 13 appears to offer a dis- 
tinct advantage to the Magnus expansion compared with 
the regular perturbation theory. On the other hand, 
evaluation of multiple integrals (9) directly is compu- 
tationally challenging, which limits the use of higher- 
order Magnus expansions for numerics. We note, how- 
ever, that the order- if universal self-refocusing condition 
C\ = . . . = Ck = is formally equivalent to 

R 1 = R 2 = . . . = R K = 0. (10) 

The matrices Rk in the latter condition are much easier to 
evaluate numerically using Eqs. (6) — (8). Yet, the ben- 
efits of the cluster theorem remain: to K-th order only 
clusters with up to K + 1 vertices need to be analyzed. 

Wc implemented the described scheme using the stan- 
dard fourth-order Rungc-Kutta algorithm for solving 
coupled differential equations, and the GSL library 15 
for matrix operations. The coefficient optimization was 
done using a combination of simulated annealing and the 
steepest descent method. The trial pulse shapes were 
encoded in terms of their Fourrier coefficients, 

V(t+r/2) = A m cos(mflt)+B m sin(mfit), (11) 

m 
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where the angular frequency Q = 2tt/t is related to the 
full pulse duration r. The target function for single- 
pulse optimization included the sum of the magnitudes 
squared of the matrix elements of the zeroth-order mis- 
match matrix [Uo(t) — ^target], and of the matrices Rk{r), 
k = 1, . . . , K. The minimization continued until these 
contributions went down to zero with the numerical pre- 
cision (typically, eight digits or more). 

As the simplest application of the formalism, we de- 
signed a number of inversion (ir-) pulse shapes 14 , self- 
refocusing to various degrees with respect to the Ising 
interaction; their coefficients are listed in TAB. II. To 
reduce the bandwidth of an arbitrary sequence of such 
pulses, we required additionally that the function (11) 
vanishes along with a number of its derivatives V^ l \t), 
I = 1, 2, . . . , 2L — 1, at the ends of the interval, t = 0, r. 

These shapes can work in known high-order pulse 
sequences 16 as a drop-in replacement of hard (or short 
Gaussian) pulses. We note that in our setup there is no 
gap between subsequent pulses, the pulses follow back to 
back with the repetition period r. The system is "fo- 
cused" at the end of each time interval. Such a scheme 
with a common "clock" time r is convenient, e.g., for par- 
allel execution of quantum gates in different parts of the 
system. For each qubit, various pulses (or intervals of no 
signal) can be executed in sequence. The performance of 
such sequences can be analyzed in the same manner as 
that of a single pulse. Namely, we integrate Eqs. (6), (7), 
and (8) over the full duration t of the pulse sequence; the 
order of the sequence is the number K of the exactly can- 
celled terms in the perturbative expansion of R(t). After 
N = t/r steps, the error in the unitary evolution matrix 
would scale as oc Nt k+1 = tr K ; the corresponding gate 
fidelity (defined as the probability of error, either average 
or maximum) would scale as 1 — 0{t 2K ). 

In Tab. I, we illustrate the quality of the obtained 
pulses by comparing their performance in several refo- 
cusing sequences for different models. "Ising" : the Ising- 
only interaction [Eq. (3) with all J 1 = 0]; "xxz": the 
xxz spin chain with both J z and J non-zero; "bath" : 
Ising spin chain coupled to a thermal bath generating 
slow (compared to r) phase modulation, simulated as 
H a [Eq. (4)] with random time-independent coefficients 
B* (see further discussion on open systems below). The 
pulse sequences are listed in the caption; these are "best" 
sequences at given length for all pulse shapes found by 
exhaustive search (high-order sequences 16 equivalent for 
hard pulses do not necessarily have equal orders here). 
The fact that such a brute-force optimization approach 
works is entirely due to the efficiency of the method. 

The most interesting is the length-8 sequence 
"X? Y$Zl Y 2 2 Y 2 2 xl Y 2 X 2 " , where X 2 is a 7r-pulse in x- 

2 

direction applied on every odd site, Y 2 is a 7r-pulse in 
negative y-direction on even sites, etc. This sequence is 
the best among the length-eight sequences for both the 
Ising and the xxz (J 1 ^ 0) models, and, additionally, it 
protects every qubit from phase decoherence due to low 
frequency noise. Our second-order self-refocusing pulses 



are clearly advantageous, especially if the Ising coupling 
is dominant. The corresponding errors scale as (J z r) 6 
compared with that for the standard (first-order) Her- 
mitian pulse where gate error scales as (J z r) 4 [the gate 
fidelities differ from unity by 0({J z t) 12 ) and C((J z r) 8 ) 
respectively] . 

These pulses were designed for use in systems with 
dominant Ising coupling, and this is the situation where 
they are most useful as a replacement of, say, Gaussian 
pulses. For example, when the pulse Qi along with anal- 
ogously designed second order ir/2 and 2tt pulses were 
used to simulate the BBi composite pulse 17 designed to 
compensate for amplitude errors to third order, the re- 
sults for a single spin were essentially identical to those 
with Gaussian pulses, with errors cubic in the amplitude 
mismatch. However, when used in an Ising chain, the 
performance of the BBi sequence with Gaussian pulses 
deteriorated linearly in J z t already with zero amplitude 
mismatch, while for our second-order pulses the addi- 
tional error was smaller, scaling as the product of (J z t) 
and the amplitude mismatch. Clearly, if the two sources 
of errors are comparable, combining high-accuracy BBi 
composite pulse and the second order pulses may be su- 
perficial; simpler pulse sequence and/or pulses with first 
order compensation could give a comparable accuracy. 

Open systems: Qualitatively, the effect of the refo- 
cusing pulses on the thermal bath coupling H [Eq. 4] can 
be most readily understood in the rotating frame defined 
by the bare evolution operator Uo(t) [Eq. 6], 

H a (t) = U^t)H a U = \HM Pl , qi )Qf {t)a£ . (12) 

For refocusing, the rotation matrices (t) are periodic 
with the full sequence period f; they can be written as a 
sum of harmonics with the main frequency fi = 2n/f, 

QT\t) = Y, c ™ e ~ iSlmt > n m ^ m n. (13) 

m 

The constant-field first-order refocusing condition (aver- 
age Hamiltonian vanishes to leading order) is equivalent 

to a cancellation of some linear combinations of C^q 

(e.g., C z n Q — for phase noise assumed for thermal bath 
in Tab. I). As a result, the environmental modes at low 
frequencies get modulated and are effectively replaced by 
those at higher frequencies, u> — > u> + f2 m , m =/= 0, leading 
to a significant reduction of the decoherence caused by 
resonant decay processes 18-20 . On the other hand, fast 
modes are mostly unaffected; modulation has essentially 
no effect on a "fast" (e.g., (5-correlated) thermal bath. 

Quantitatively, the effect of refocusing can be under- 
stood with the quantum kinetic equation (QKE) in the 
rotating frame, with the kernel accurate at least to or- 
der K to analyze order- if refocusing 21 . For large 17, 
the density-matrix dynamics separates onto sectors with 
frequencies around Cl m . The slow sector, m — 0, car- 
ries the main part of the total weight, with that of the 
remaining (generally, rapidly-decaying) sectors totalling 
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Ai 


A 2 


A 3 


Aa 


As 


Si 


0.5 


-1.2053194466 


0.4796460175 


0.2256734291 






s 2 


0.5 


-1.1950755990 


0.7841246569 


0.0738054432 


-0.1628545011 




Qi 


0.5 


-1.1374003264 


1.5774784244 


-0.6825954606 


-0.2574826374 




Q2 


0.5 


-1.0965122417 


1.5309957409 


-1.1470791601 


0.0020722004 


0.2105234605 



TABLE II: Fourrier coefficients for the constructed pulses, see Eq. (11). Shapes Sl and Ql respectively are first (K = 1) and 
second (K = 2) order self-refocusing inversion pulses for the Ising spin coupling; the fixed-time errors scale with the duration 
of the pulse as oc (tJ z ) k [cf. Eq. (3)]. These shapes have 2L derivatives vanishing at the ends of the interval. 



- A(0)/O 2 , where A(t-t') = \\(B^(t)B^' (t'))\\ is a norm 
of the correlation matrix of the fluctuating field. Only 
the dynamics in the slow sector is protected by the re- 
focusing. In particular, the analysis of the QKE with 
the leading second-order kernel shows that already with 
first-order (K = 1) constant-field refocusing direct de- 
cay processes require excitations at frequencies w > f2, 
which may dramatically reduce the dissipative part of 
the QKE kernel. The non-resonant reactive processes 
are also suppressed: the rate of phase errors is ~ A(0)/f2 
with K = 1 refocusing and ~ |A"(0)|/O 3 (primes denote 
time derivatives) with K = 2 refocusing, as, e.g., for 
length- 8 sequence in Tab. I. Generally, these results 21 
apply equally for soft- and hard-pulse refocusing, and 
are consistent with established results on kinetics of few- 
level systems in r.f. field 18 , and with the properties of 



hard pulse sequences for low- frequency environment 20,22 . 

To conclude, we presented an efficient scheme for de- 
signing high order soft pulses and soft-pulse sequences in 
a scalable fashion, without the need for solving the full 
Hamiltonian. Soft (narrow-band) pulses are indispens- 
able for their selectivity and reduced coupling to envi- 
ronmental modes, which in turn suppresses signal distor- 
tions and heating. Use of high order pulses is especially 
efficient if one interaction (e.g., the Ising term) is dom- 
inant. High order pulse sequences generally offer better 
accuracy and can dramatically reduce the decoherence 
due to coupling with low- frequency environment. 
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